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Abstract 

It is shown that the critical layer analysis, involved in the linear theory of internal 
modes, can be extended continuously into the early nonlinear regime. For the m = 1 
resistive mode, the dynamical analysis involves two small parameters: the inverse of 
the magnetic Reynolds number S and the m = 1 mode amplitude A, that measures 
the amount of nonlinearities in the system. The location of the instantaneous critical 
layer and the dominant dynamical equations inside it are evaluated self-consistently, 
as A increases and crosses some 5-dependent thresholds. A special emphasis is put 
on the influence of the initial (/-profile on the early nonlinear behavior. Predictions 
are given for a family of ^-profiles, including the important low shear case, and 
shown to be consistent with recent experimental observations. 
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The m = n = 1 internal modes, such that the safety factor goes below one for 
some inner radius, remain critical macroscopic modes for large scale tokamak 
plasma dynamics and confinement. They are particularly involved in sawtooth 
oscillations and crashes. This is a common deleterious phenomenon as conven- 
tional tokamak discharges eventually operate with qo < 1 since current density 
tends to a peaked profile. Additionally, the m — n — 1 internal modes form 
a laboratory prototype for reconnection. Such phenomena typically proceed 
beyond linear regime. 

We shall consider here the m — n — 1 purely resistive mode [1] that is ideally 
marginally stable. The original motivation of this work was to understand the 
growth of the m = 1 resistive mode up to its nonlinear saturation, on the basis 
of some striking numerical simulations performed by Aydemir [2] and previous 
observations [3]. Within the reduced MHD framework in cylindrical coordi- 
nates and some given g-profile [2], the time behavior of the kinetic energy in 
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the m = 1 mode amounts to an initial exponential growth consistent with the 
linear regime, followed by a transient stage where the growth rate decreases, 
that is brutally interrupted by a second exponential growth in the nonlinear 
regime. This second exponential stage eventually terminates, as the kinetic 
energy in the m — 1 mode saturates which coincides with the completion of 
magnetic reconnection. 

The reduced MHD system under consideration reads 



Helical symmetry is assumed: The poloidal and toroidal angles, respectively 6 
and ip, only come in through the helical angle a = ip—Q. and ip are the plasma 
velocity and helical magnetic field potentials: the velocity is v = <p x Vj_0 and 
the magnetic field B = B Oip + ip x Vj_ {ip — r 2 /2). U = V^_0 is the vorticity 
and J = V 2 ± ip the helical current density, with = r~ x d r rd r + r~ 2 d 2 a . 
Poisson brackets are defined by [0, U\ = —(p ■ (Vj_0 x V±U) = r^ 1 (d r (f)d a U — 
d r Ud a 4>). Eqs. (l)-(2) are dimensionless: Time has been normalized to the 
poloidal Alfven time, the radial variable r to the minor radius, and 77 is the 
inverse of the magnetic Reynolds number S, and is given by the ratio of the 
poloidal Alven time to the resistive one. In high-temperature fusion plasmas, 
i] is typically much smaller than one. 

Consider equilibria such that, for some internal radius r s o < 1, q{r s o) = 1, 
that is ip' (r s o) = 0. Then, due to the Ohm's law (2), plasma volume divides 
in two region. Far from the q — 1 surface (outer domain), plasma behaves 
ideally whereas, in the vicinity of r s0 (inner region), resistivity plays a crucial, 
destabilizing, role. Linear theory [1] uses asymptotic matching analysis to 
provide m — 1 eigenfunctions in the form A(t)fi,(r) exp(ia) valid in the whole 
domain. In the outer (ideal) domain, this solution is valid, that is nonlinear 
effects are negligible, as long as i < 1 [4]. Injecting the linear solutions 
ipi(r,a,t) = A{t)ip L {r-) expect) and 0i(r, a,t) = A(t)0 £ (r) exp(ia) into (1)- 
(2) calls for an amplitude expansion. The procedure has been given in Refs. 
[4,5]. The particularity of the linear radial eigenfunctions ipL{ r ) and 4>L{ r ), 
that needs a careful consideration, is that they have strong gradients inside 
the critical layer. More precisely, their radial derivatives are of the order of the 
inverse of the critical layer width, that is O^ 1 ^). This means in particular 
that this approach restricts to situations strictly above marginal stability and 
where the linear regime is well defined, with clear scalings, yielding the resistive 
ordering, and non-pathological g-profiles (in the sense of Ref . [6]). 

We wish then to answer the question: "How does the m — 1 resistive mode 
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develop into the nonlinear regime ?" To do this, let us first recognize that the 
problem involves two small parameters. An obvious one is the resistivity rj. 
However, considering it to be the only one small parameter, in some perturba- 
tion analysis with conventional expansions of the type / = fo+rjfi + . . . would 
lead to a dead end: this would bring up a singular expansion, with additional 
77 In (77) terms, with no asymptotic validity unless assuming that the mode am- 
plitude is always kept vanishingly small. It is interesting to note that such a 
procedure would actually be valid for the tearing mode with the small param- 
eter limit A' [7,8]. In the present case, such a perturbation analysis would be 
ill-posed. A second small parameter enters the game, the m = 1 mode ampli- 
tude A that can be viewed as an indicator of the amount of nonlinearities in 
the system. As previously said, the approach will then be that of an amplitude 
expansion. 

The first step will be to determine the end of validity of the linear regime. In 
the outer domain, this occurs for A of order one [4] but, in the inner domain, 
the linear solution breaks earlier. This occurs when mode coupling terms such 
as [0i, Ui] becomes of the same order order as linear terms, that is for A > rf^ . 
At this point, m = and m = 2 components begin to be "fed" nonlinearly 
by mode coupling terms: the m = and m = 2 modes are nonlinearly driven. 
However, these mode coupling terms, quadratic in A, do not affect the m — 1 
dynamics so that one could say that the m = 1 mode is still linear. At this 
stage, it is easy to check that the dominant equations on the m — 1 component 
are still the linear ones. This means that the radial structure of the solution 
should remain close to the linear one. Given that, it is possible to include the 
correction to the linear theory due to the new location of the critical layer. 
Because of the m — 1 perturbation, the critical layer is not expected to remain 
fixed at r s0 . The real helical magnetic field potential inside the critical layer 
(for \r — r s0 \ < rf^) is 



In writing down the critical layer equations, the instantaneous surface r s (a, t), 
defined by d r ip [r s (a,t)} = 0, is important as the location where dynamical 
equations turn singular and non-ideal effects come into play. As mode cou- 
plings do not affect the second order m — 1 dynamics, one can keep the linear 
m = 1 radial structure but introduce the corrections due to the motion of 
the d r ip = surface. This yields a differential equation [4,5] for the m — 1 
amplitude A(t) valid below the onset of "truly nonlinear" cubic nonlinearities 
on m — 1. This is given by 
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with 



7(*) 7(* = 0) 1l 

'frMl\ / ^[r 3 (a,t=0)} \ <(r s(1 ) ' 

> r 3 (a,t) I a \ r s (a,t=0) / a r aQ 



(5) 



where (-) a = {2tt)~ 1 J" 27r -da denotes the m = average. As shown in Refs. 
[4,5], the differential equation (4), with ^(t) given by (5), may actually be 
approximated by the quadratic expression 

^ = lL A(t)-C A 2 (t), (6) 



with 



Co 



q'(r s0 ) + r s0 q"(r s0 ) - 2r s0 q'(r s0 f 
V2nr 2 s0 q'(r s0 ) 



(7) 



It may be useful to remind here that the problem has been rendered dimen- 
sionless. The solution of (6) is 

A(t) = A ° eXP {lht) (8) 

U 1 + A C / lL [exp ( lL t) - 1] { ) 



that tends to 7l/Co as t — > oo. It is interesting to note that the curve A/A has 
a universal form depending only on the rescaled time 7^t and on the parameter 
AqCq/^l containing the g-profile properties at r s0 . It is also interesting to note 
that, at time U = 7^ In [7^/ (A Cq) — 1], A(t) possesses an inflexion point so 
that, around that time, the effective behavior of A is approximately algebraic 
(being linear). All this assumes that Cq is positive. A negative Co would yield 
a transient explosive faster-than-exponential behavior. However, the validity 
of (6) is limited because some cubic nonlinearities should come into play. 
For S not too large, so that the instantaneous second order location of the 
critical layer has some large overlap with the initial linear critical layer at r s0 
around the X-point, it is possible [4,5] that those cubic terms show up in a 
spectacular manner. When collecting terms cubic in A, it turns out that in 
this overlap domain the convective derivative due to the motion of the critical 
layer dominates the ordinary time derivative and equilibrates mode coupling 
terms yielding, together with (6), an effective amplitude equation of the form 

with c = Oil) some constant and 7^ the growth rate reached by the m = 1 
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Fig. 1. Different (/-profiles of the family q(r) = qo |l + 
go = 0.9 and g a = 3 for A = 0.6 (dashed line), A = 2 (plain), A 
A = 30 (dot-dashed line). The bold line is the q = 1 threshold. 
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mode at the onset of cubic terms. This takes place for A 
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Let us finally investigate the influence of the g-profile in the onset of the nonlin- 
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ear regime. The g-pronles q(r) = qo jl + r 
by A, can depict different possible experimental situations (see Figure 1). The 
lowA case is consistent with a very peaked current profile. Such a kind of pro- 
file was used for instance by Biskamp in his 1991's simulations [9] of the same 
system. He observed a transition from the linear exponential m = 1 growth 
towards an algebraic behavior. On the contrary, the large A case coincides with 
a low shear situation with a flat current profile within the q — 1 radius. This is 
reminiscent of recent experimental investigations undertaken e.g. in JET un- 
der the "hybrid" scenario, with a wide area of low magnetic shear and central 
safety factor close to and below one. Buratti and coworkers have reported in 
this case the wide emergence of "slow sawteeth" [10,11] where the mode has 
the same spatial structure as the kink-like sawtooth precursor with n — 1 but 
grows very slowly and enters the nonlinear regime at the linear growth rate. 
Although this case may not strictly correspond to the purely resistive mode, 
the present analysis on the onset of nonlinear effects should be transposable. 
Figure 2 may indeed propose an explanation for these observations. In the 
case of a very peaked g-profile (e.g. with A = 0.6), the integration of m = 1 
amplitude evolution Eq. (8) before the onset of third order convective effects 
shows that A should saturate before this third order threshold A ~ r\ l l 2 is 
reached. In particular, this may explain Biskamp's observations [9] of a tran- 
sition to an algebraic stage with no subsequent nonlinear exponential growth. 
The A = 2 case just corresponds to the g-profile taken by Aydemir in Ref. [2]. 
Here the third order convective stage showing a second stage of exponential 
growth can be reached. It is clear from Fig. 2 that the nonlinear growth rate 
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Fig. 2. Evolution (in linear-log scale) of the amplitude of the m = 1 mode, given 
by Eq. (8), as a function of time normalized to the linear growth rate for the 
^-profiles displayed in Fig. 1 (with the same plot styles) before the onset of "third 
order" nonlinear regime. The bold horizontal line marks the threshold of third order 
convective terms for A ~ S~ l l 2 for the value S = 10 7 . The initial amplitude is 
A Q = 1.2 x 1(T 5 . 

7at/,, that is the growth rate of the m — 1 mode when A crosses the threshold 
(A ~ if-l 2 ) is (slightly) smaller than 7^ as A has turned bending. Indeed the 
simulations of Ref. [2] give the numerical value 7tvl — 7l/2. Yet, for larger 
values of A, corresponding to a g-profile of the kind studied by Buratti et al. 
[10,11], it is clear from Figure 2 that ^ NL would be almost equal to 7 L : This 
is an explanation of the fact that, in "slow sawteeth", the m — 1 mode enters 
the nonlinear regime at almost the linear growth rate. 



Acknowledgements 



The author is greatly indebted to B. Coppi for his advice and explanations 
and thanks the organizers and participants of the 9th Plasma Easter Meeting 
held in Torino for the nice meeting. 



References 



[1] B. Coppi, R. Galvao, M.N. Rosenbluth and P.H. Rutherford, Soviet Journal of 
Plasma Physics 2, 3276 (1976). 

[2] A.Y. Aydemir, Phys. Rev. Lett. 78, 4407 (1997). 

[3] B.V. Waddell, M.N. Rosenbluth, D.A. Monticello and R.B. White, Nucl. Fusion 
16, 3 (1976). 



6 



[4] M.C. Firpo, Phys. Plasmas 11, 970 (2004). 

[5] M.C. Firpo and B. Coppi, Phys. Rev. Lett. 90, 095003 (2003). 

[6] R. Fitzpatrick, Plasma Phys. Control. Fusion 31, 1127 (1989). 

[7] D.F. Escande and M. Ottaviani, Physics Letters A 323, 278 (2004). 

[8] F. Militello and F. Porcelli, Physics of Plasmas 11, L13 (2004). 

[9] D. Biskamp, Phys. Fluids B 3, 3353 (1991). 

[10] P. Buratti, B. Alper, A. Becoulet, P. Belo, C. Gormezano, P. Smeulders and 
EFDA-JET Contributors, 31st EPS Conference on Plasma Physics, London, 
Europhysics Conference Abstracts 28G, P. 1-165 (2004). 

[11] C. Gormezano, A. Becoulet, P. Buratti et al., Plasma Physics and Controlled 
Fusion 46, B435 (2004). 



7 



